Constant-time solution to the global optimization problem using Brüschweiler’s ensemble search algorithm

A constant-time solution of the continuous global optimization problem (GOP) is obtained by using an ensemble algorithm. We show that under certain assumptions, the solution can be guaranteed by mapping the GOP onto a discrete unsorted search problem, whereupon Brüschweiler’s ensemble search algorithm is applied. For adequate sensitivities of the measurement technique, the query complexity of the ensemble search algorithm depends linearly on the size of the function’s domain. Advantages and limitations of an eventual NMR implementation are discussed. PACS numbers: 03.67.Lx, 33.25.+k, 76.60.−k 1. The global optimization problem Optimization problems are ubiquitous and extremely consequential. Their theoretical and practical interest has continued to grow from the first recorded instance of Queen Dido’s problem [1] to present day forays into complexity theory and large-scale logistics applications (see [2–6] and references therein). The formulation of the global optimization problem (GOP) is deceptively simple: find the absolute minimum of a given function—called the objective function—over the allowed range of its variables. Sometimes, the function whose global minimum is to be found is not specified in analytic form and must be evaluated by a computer program. When we can only access the output of the computation, the program acts as a black-box tool, called an oracle. Of course, the brute force approach of evaluating the function on its whole domain is either impossible if the variables are continuous, or prohibitively expensive if the variables are discrete, but take values in large domains and in high-dimensional spaces. Since, in general, each oracle invocation (function evaluation) involves an expensive computational process, the 0305-4470/03/240399+09$30.00 © 2003 IOP Publishing Ltd Printed in the UK L399 L400 Letter to the Editor number of function evaluations needs to be kept to a minimum. It is not surprising that, together with accuracy, this number probably provides the paramount criterion in comparing the efficiency of competing optimization algorithms. The primary difficulty in solving GOPs stems from the fact that the condition for determining minima, namely annulment of the gradient of the objective function, is only necessary (the function may have another type of critical point) and local, i.e. it does not distinguish between local and global minima. The generic strategy to find the global minimum involves two main operations, namely: (i) descent to a local minimum and (ii) search for a new descent region. However, this generic strategy is powerless for certain problems. The following one-dimensional example (the golf course problem) illustrates the difficulties. Define the function f : [0, 1] by f (x) = { 0 for 0 x a − /2 and a + /2 x 1 −1 for a − /2 x a + /2 (1) where 0 < 1, and a is a point in the interval ( /2, 1 − /2), but is otherwise unknown. From this definition, it is clear that to find the absolute minimum of f , one should evaluate it within the well-interval around the unknown point a. If this function is defined like an oracle, i.e. if one does not know its definition above and the value of the number a, the probability P of choosing the variable x within the well-like interval where f (x) = −1, is . In the n-dimensional version of this problem, the probability P becomes , so the complexity of the problem grows exponentially with n (the dimensionality curse). Moreover, knowledge about the (partial) derivative(s) of f would not help, since they are all zero whenever defined. In the light of this example, it seems that without additional information about the structure of the function there is no reasonable expectation of deciding upon an efficient optimization strategy and one is left with either strategies that have limited applicability or the exhaustive search option. Thus, new approaches are needed, which take advantage of additional information to reduce the complexity of the problem to a manageable level. We emphasize that this information is actually available in some classes of GOP—for instance it happens to be available for the golf course problem above—but cannot be taken advantage of within conventional optimization algorithms. Here we present an approach that uses this information efficiently to map the continuous GOP into a discrete search problem (DSP). Once this mapping is completed, we apply a variant of Brüschweiler’s search algorithm, as described in the following section, to obtain a value within the basin of attraction of the global minimum in a logarithmic number of function evaluations that scales logarithmically with respect to the cardinality of the resulting discrete set. While the parallelism of ensemble computing provides an ‘exponential speed-up’, it also demands an exponentially large number of processors. Thus, we use the phrase ‘constant-time’ to account for the relationship between speed and resources. 2. Mapping the continuous GOP to a discrete search problem Consider a real function of d variables, f (x), x = (x1, x2, . . . , xd) that has a certain degree of smoothness. The precise definition is not critical for the argument below. Without restricting the generality, we can assume that f is defined on [0, 1] → [0, 1]. Assume now that: (i) there is a unique global minimum and its value is zero; (ii) there are no local minima whose value is infinitesimally close to zero; in other words, the values of the other minima are larger than a constant δ > 0 and (iii) the size of the basin of attraction for the global minimum measured at height δ is known to be rδ . We note that these assumptions can be actually verified for Letter to the Editor L401 the problem at hand, if additional information is available. Also, these assumptions may in general be relaxed to some extent, as outlined in the first part of section 4. The implementation idea is the following. Instead of f (x), consider the function g(x) := (f (x))1/m. For sufficiently large m, this function will take values very close to 1, except in the vicinity of the global minimum, which will maintain its original value, namely zero. Of course, other transformations may be used to achieve essentially the same result. Now divide the hypercube [0, 1] in small d-dimensional hypercubes with sides x = 1/M , where M is a positive integer. At the midpoint of each of these hypercubes, define the function h(x) := 1 − INT [g(x) + 1/2] (INT denotes the integer part). As a result of the transformations above, we reduced the GOP to a DSP for h : {1, 2, . . . , N} → {0, 1}, which is an integer-valued function defined on a discrete set of N := 2n+1 points. The function h(i) is known to be zero for all inputs i, except for one special input i = q , where h(q) = 1. Here N = M and the domain on which the function h is zero depends on m and on the size, rδ , of the basin of attraction of the global minimum at height δ. Thus the problem becomes to find the value of the special input q efficiently. Recently, we proposed an adaptation of Grover’s quantum search algorithm to solve the continuous GOP [7], which requires O( √ N) function evaluations [8–10]. In this letter, we propose to further improve the efficiency, by using an ensemble search algorithm, which is described in detail in the following section. The algorithm is adapted from Brüschweiler’s ensemble search algorithm [11, 12], using a parallel technique to achieve exponential speedups for problems up to a certain size, dictated by the sensitivity of the measurement process. As mentioned before, efficiency is understood in relation to the query complexity of the algorithm. Indeed, when N is large and the function evaluation (i.e. the oracle) is costly in terms of computational complexity, reducing the number of function evaluations is critical. Application of the ensemble search algorithm to the function h results in a point that returns the value zero. It is easy to see that, by construction, this point belongs to the basin of attraction of the global minimum. We return then to the original function f and apply the descent technique of choice that will lead to the global minimum. If the basin of attraction of the global minimum is narrow, the gradients of the function f may reach very large values which may cause overshots. Once that phase of the algorithm is reached, one can apply a scaling (dilation) transformation that maintains the descent mode but moderates the gradients. On the other hand, as one approaches the global minimum, the gradients become very small and certain acceleration techniques based on non-Lipschitzian dynamics may be required [18, 19]. If the global minimum is attained at the boundary of the domain, the algorithm above will find it without additional complications. 3. Brüschweiler’s ensemble search algorithm Following Grover and Shor’s breakthrough algorithms for quantum computing [8, 13, 14], an alternative paradigm for computing has been suggested by Madi, Brüschweiler and Ernst, which operates on ensembles, i.e., mixed states of identical spin sytems, using a spin Liouville space formalism for density operators [11]. The Liouville space for a system of n spins contains the n × n density matrices representing all the possible quantum states of the spins. Ensemble algorithms are not quantum algorithms, insofar as they do not involve entanglement of quantum states, so they only use a small subset of the Liouville space. Throughout this letter, ‘mixed states’ describe a statistical ensemble of spins, rather than individual spins, so that each element of the ensemble performs part of the computation, in the same way as each processor in a classical parallel computer. L402 Letter to the Editor This new paradigm exploits the parallelism offered by simultaneously acting on many different input states in an ensemble of spins. In contrast, quantum computing with pure states relies on the parallelism of entangled states, to perform operations in a Hilbert space which is the tensor product of multiple qubits. The advantage of ensemble algorithms is that they may be exponentially faster than the equivalent quantum algorithms, for adequate measurement sensitivities. Two realizations of ensemble algorithms are provided by Brüschweiler’s search algorithm [12] and the summing algorithm proposed by the authors [15], which illustrate the trade-off between memory and speed capabilities. Brüschweiler’s algorithm for searching an unsorted database uses the ensemble computing paradigm in the context of NMR technology. The strategy employs binary partition of the N elements in the database, to find the desired element after O(log2 N) oracle queries. We note that the lower bound derived by Nayak and Wu [9] for the query complexity of quantum search algorithms does not apply here, since the quantum lower bound [10] used in their derivation assumes that the initial state is a pure state. Ensemble computing requires an exponentially larger set of physical resources to encode the same number of distinct input states compared to quantum computing with pure states. On the other hand, ensemble computing holds the important advantage that it is insensitive to the decoherence time of the spins, which is an outstanding limiting factor for quantum computations involving entangled states. A ‘divide and conquer’ scheme is used to test whether the special input q belongs to exponentially finer and finer partitions of the set of input values. The algorithm is envisaged for a physically realizable system consisting of binary-valued spins, therefore the number of spins needed in the input register varies from n to 1, to accommodate the number of input values in partitions ranging in size from N/2 values down to 2 values, respectively. It is also assumed that a 1-spin output register is available, to encode the value of h ∈ {0, 1}. The first step of the algorithm is to divide the input values into two equal partitions, χ− 1 = {1, . . . , N/2} and χ+ 1 = {N/2 + 1, . . . , N}, and test whether the special input q belongs to χ− 1 or χ + 1 . The numerical subscript, k = 1, indicates the size of the partition as a fraction (1/2) of the entire set of input values, χ = {1, . . . , N}, and the superscript ± differentiates between the lower and upper partitions. The n-spin input register is initialized as an ensemble, namely an equally-weighted mixed state,